
clear
rng default

%% Empirical distribution from data on cohort 1950-1951

%%Men%%
year=1950;
% Probability distribution of Ychosen|X
condprobsmen_temp=load('probs_allyears_data.mat',['condprobsmen' num2str(year) '_140']);
condprobsmen=condprobsmen_temp.(['condprobsmen' num2str(year) '_140']);
PYX_cond=[condprobsmen(:,2:end) condprobsmen(:,1)]; %5x6
%Order: P(Ychosen=1|X=1) P(Ychosen=2|X=1) P(Ychosen=3|X=1) P(Ychosen=4|X=1) P(Ychosen=5|X=1) P(Ychosen=0|X=1) 
%       P(Ychosen=1|X=2) P(Ychosen=2|X=2) P(Ychosen=3|X=2) P(Ychosen=4|X=2) P(Ychosen=5|X=2) P(Ychosen=0|X=2) 
%       P(Ychosen=1|X=3) P(Ychosen=2|X=3) P(Ychosen=3|X=3) P(Ychosen=4|X=3) P(Ychosen=5|X=3) P(Ychosen=0|X=3) 
%       P(Ychosen=1|X=4) P(Ychosen=2|X=4) P(Ychosen=3|X=4) P(Ychosen=4|X=4) P(Ychosen=5|X=4) P(Ychosen=0|X=4) 
%       P(Ychosen=1|X=5) P(Ychosen=2|X=5) P(Ychosen=3|X=5) P(Ychosen=4|X=5) P(Ychosen=5|X=5) P(Ychosen=0|X=5) 
probsmen_temp=load('probs_allyears_data.mat',['probsmen' num2str(year) '_140']);
PX=probsmen_temp.(['probsmen' num2str(year) '_140']); %5x1
%Order: P(X=1) P(X=2) P(X=3) P(X=4) P(X=5) 

%Regroup to have 2 types plus singles
PYX_joint=PYX_cond.*kron(PX, ones(1,6));
PYX_joint_new_1=[PYX_joint(:,1)+PYX_joint(:,2)  PYX_joint(:,3)+PYX_joint(:,4)+PYX_joint(:,5)  PYX_joint(:,6)];
PYX_joint_new=[PYX_joint_new_1(1,:)+PYX_joint_new_1(2,:); PYX_joint_new_1(3,:)+PYX_joint_new_1(4,:)+PYX_joint_new_1(5,:)]; 
PX_emp=[PX(1)+PX(2); PX(3)+PX(4)+PX(5)];
PYX_cond_emp=PYX_joint_new./(kron(PX_emp, ones(1,3)));
PYX_cond_emp=[PYX_cond_emp(1,:) PYX_cond_emp(2,:)];


%%Women%%
year=1951;
% Probability distribution of Xchosen|Y
condprobswomen_temp=load('probs_allyears_data.mat',['condprobswomen' num2str(year) '_140']);
condprobswomen=condprobswomen_temp.(['condprobswomen' num2str(year) '_140']);
PXY_cond=[condprobswomen(:,2:end) condprobswomen(:,1)]; %5x6
%Order: P(Xchosen=1|Y=1) P(Xchosen=2|Y=1) P(Xchosen=3|Y=1) P(Xchosen=4|Y=1) P(Xchosen=5|Y=1) P(Xchosen=0|Y=1) 
%       P(Xchosen=1|Y=2) P(Xchosen=2|Y=2) P(Xchosen=3|Y=2) P(Xchosen=4|Y=2) P(Xchosen=5|Y=2) P(Xchosen=0|Y=2) 
%       P(Xchosen=1|Y=3) P(Xchosen=2|Y=3) P(Xchosen=3|Y=3) P(Xchosen=4|Y=3) P(Xchosen=5|Y=3) P(Xchosen=0|Y=3) 
%       P(Xchosen=1|Y=4) P(Xchosen=2|Y=4) P(Xchosen=3|Y=4) P(Xchosen=4|Y=4) P(Xchosen=5|Y=4) P(Xchosen=0|Y=4) 
%       P(Xchosen=1|Y=5) P(Xchosen=2|Y=5) P(Xchosen=3|Y=5) P(Xchosen=4|Y=5) P(Xchosen=5|Y=5) P(Xchosen=0|Y=5) 
probswomen_temp=load('probs_allyears_data.mat',['probswomen' num2str(year) '_140']);
PY=probswomen_temp.(['probswomen' num2str(year) '_140']); %5x1
%Order: P(Y=1) P(Y=2) P(Y=3) P(Y=4) P(Y=5) 

%Regroup to have 2 types plus singles
PXY_joint=PXY_cond.*kron(PY, ones(1,6));
PXY_joint_new_1=[PXY_joint(:,1)+PXY_joint(:,2)  PXY_joint(:,3)+PXY_joint(:,4)+PXY_joint(:,5)  PXY_joint(:,6)];
PXY_joint_new=[PXY_joint_new_1(1,:)+PXY_joint_new_1(2,:); PXY_joint_new_1(3,:)+PXY_joint_new_1(4,:)+PXY_joint_new_1(5,:)]; 
PY_emp=[PY(1)+PY(2); PY(3)+PY(4)+PY(5)];
PXY_cond_emp=PXY_joint_new./(kron(PY_emp, ones(1,3)));
PXY_cond_emp=[PXY_cond_emp(1,:) PXY_cond_emp(2,:)];

%% Get the systematic match surplus under the Logit assumption
%Notation: u12=man is type 2 and chooses woman of type 1
%u11=u1(1), u12=u1(2), u22=u2(2), u21=u2(1);
u11_emp=log(PYX_cond_emp(1)/PYX_cond_emp(3)); %u1(1)=log(P(Ychosen=1|X=1)/P(Ychosen=0|X=1))
u12_emp=log(PYX_cond_emp(4)/PYX_cond_emp(6)); %u1(2)=log(P(Ychosen=1|X=2)/P(Ychosen=0|X=2))
u22_emp=log(PYX_cond_emp(5)/PYX_cond_emp(6)); %u2(2)=log(P(Ychosen=2|X=2)/P(Ychosen=0|X=2))
u21_emp=log(PYX_cond_emp(2)/PYX_cond_emp(3)); %u2(1)=log(P(Ychosen=2|X=1)/P(Ychosen=0|X=1))
%Notation: v12=woman is type 2 and chooses man of type 1
%v11=v1(1), v12=v1(2), v22=v2(2), v21=v2(1); 
v11_emp=log(PXY_cond_emp(1)/PXY_cond_emp(3)); %v1(1)=log(P(Xchosen=1|Y=1)/P(Xchosen=0|Y=1))
v12_emp=log(PXY_cond_emp(4)/PXY_cond_emp(6)); %v1(2)=log(P(Xchosen=1|Y=2)/P(Xchosen=0|Y=2))
v22_emp=log(PXY_cond_emp(5)/PXY_cond_emp(6)); %v2(2)=log(P(Xchosen=2|Y=2)/P(Xchosen=0|Y=2))
v21_emp=log(PXY_cond_emp(2)/PXY_cond_emp(3)); %v2(1)=log(P(Xchosen=2|Y=1)/P(Xchosen=0|Y=1))          

Phi11_emp=u11_emp+v11_emp;
Phi12_emp=u21_emp+v12_emp;
Phi21_emp=u12_emp+v21_emp;
Phi22_emp=u22_emp+v22_emp;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Primitives
nm=6000; %number of men
nw=6000; %number of women
supX=[1 2]; %support of X (men's types)
supY=[1 2]; %support of Y (women's types)
supXY=sup(supX, supY); %support of (X,Y)
%supXY order:
%(1,1)
%(2,1)
%(1,2)
%(2,2)
PhiXY=[Phi11_emp Phi12_emp; Phi21_emp Phi22_emp];

ntypes_m=size(supX,2);
ntypes_w=size(supY,2);

%% Generate data
%Men
X=discrete(nm, round(PX_emp,3), supX); %nmx1; 
epsilon = -evrnd(0,1,nm,ntypes_w+1); %nmxntypes_w+1; 
%epsilon=(epsilon1, epsilon2, ..., epsilon0) i.i.d. Gumbel with location 0 and scale 1, independent of X
%Gumbel defined as in Wikipedia (Matlab flips the signs)

%Women
Y=discrete(nw, round(PY_emp,3), supY); %nwx1; 
eta = -evrnd(0,1,nw,ntypes_m+1); %nwxntypes_m+1; 
%eta=(eta1, eta2, ..., eta0) i.i.d. Gumbel with location 0 and scale 1, independent of Y

%% Equilibrium
%(1) Compute total surplus for each possible couple
Phi=zeros(nm,nw); %nmxnw
for j=1:nw %for each woman
    for i=1:nm %for each man
        Phi(i,j)=PhiXY(X(i), Y(j))+epsilon(i,Y(j))+eta(j,X(i));
    end
end

%(2) Reshape Phi into a vector (nm*nw)x1 following the order
%man1 woman1
%man2 woman1
%...
%man1 woman2
%man2 woman2
%...
%man1 woman3
%man2 woman3
%...
Phir=reshape(Phi,nm*nw,1); %(nm*nw)x1

%(3) Augment Phir for singles following the order
%man1 woman1
%man1 woman2
%man1 woman3
%...
%man2 woman1
%man2 woman2
%man2 woman3
%...
%man1 single
%man2 single
%...
%woman1 single
%woman2 single
%woman3 single
%...

Phira=[Phir; epsilon(:,end); eta(:,end)]; %(nm*nw+nm+nw)x1

%(3) Find equilibrium matching 
mu=equilibrium(nw, nm,Phira); %(nm*nw+nm+nw)x1


%% Conditional probabilities
% Men
Ychosen=zeros(nm,1); %nmx1 listing the type of woman man i is matched with
for i=1:nm
    %select from mu the (nw+1)x1 sub-vector reporting the choice of man i
    index=(1:nm:(nm*nw+1)).'+(i-1);
    submu=mu(index);
    %save the type of woman man i is matched with
    if submu(end)==0 %man i has not chosen to be single
        Ychosen(i)=Y(logical(submu(1:end-1)));
    else %man i has chosen to be single
        Ychosen(i)=0;
    end
end

% Probability distribution of Ychosen|X
PX=[sum(X==1)/nm sum(X==2)/nm];
PYX=[sum(X==1 & Ychosen==1)/nm sum(X==1 & Ychosen==2)/nm sum(X==1 & Ychosen==0)/nm ...
     sum(X==2 & Ychosen==1)/nm sum(X==2 & Ychosen==2)/nm sum(X==2 & Ychosen==0)/nm];
PYX_cond=PYX./kron(PX,ones(1,size(supY,2)+1)); 
%order: P(Ychosen=1|X=1) P(Ychosen=2|X=1) P(Ychosen=0|X=1) P(Ychosen=1|X=2) P(Ychosen=2|X=2) P(Ychosen=0|X=2)


% Women
Xchosen=zeros(nw,1); %nwx1 listing the type of man woman j is matched with
for j=1:nw
    %select from mu the (nm+1)x1 sub-vector reporting the choice of woman j
    index=[(1:nm).'+(j-1)*nm; nm*nw+nm + j];
    submu=mu(index);
    %save the type of man woman j is matched with
    if submu(end)==0 %woman j has not chosen to be single
        Xchosen(j)=X(logical(submu(1:end-1)));
    else %woman j has chosen to be single
        Xchosen(j)=0;
    end
end

% Probability distribution of Xchosen|Y
PY=[sum(Y==1)/nw sum(Y==2)/nw];
PXY=[sum(Y==1 & Xchosen==1)/nw sum(Y==1 & Xchosen==2)/nw sum(Y==1 & Xchosen==0)/nw ...
     sum(Y==2 & Xchosen==1)/nw sum(Y==2 & Xchosen==2)/nw sum(Y==2 & Xchosen==0)/nw];
PXY_cond=PXY./kron(PY,ones(1,size(supX,2)+1)); 
%order: P(Xchosen=1|Y=1) P(Xchosen=2|Y=1) P(Xchosen=0|Y=1) P(Xchosen=1|Y=2) P(Xchosen=2|Y=2) P(Xchosen=0|Y=2)

save('PYX_cond.mat', 'PYX_cond')
save('PXY_cond.mat', 'PXY_cond')

%% Systematic payoff  
%Notation: u12=man is type 2 and chooses woman of type 1
%u11=u1(1), u12=u1(2), u22=u2(2), u21=u2(1);
u11=log(PYX_cond(1)/PYX_cond(3)); %u1(1)=log(P(Ychosen=1|X=1)/P(Ychosen=0|X=1))
u12=log(PYX_cond(4)/PYX_cond(6)); %u1(2)=log(P(Ychosen=1|X=2)/P(Ychosen=0|X=2))
u22=log(PYX_cond(5)/PYX_cond(6)); %u2(2)=log(P(Ychosen=2|X=2)/P(Ychosen=0|X=2))
u21=log(PYX_cond(2)/PYX_cond(3)); %u2(1)=log(P(Ychosen=2|X=1)/P(Ychosen=0|X=1))
%Notation: v12=woman is type 2 and chooses man of type 1
%v11=v1(1), v12=v1(2), v22=v2(2), v21=v2(1); 
v11=log(PXY_cond(1)/PXY_cond(3)); %v1(1)=log(P(Xchosen=1|Y=1)/P(Xchosen=0|Y=1))
v12=log(PXY_cond(4)/PXY_cond(6)); %v1(2)=log(P(Xchosen=1|Y=2)/P(Xchosen=0|Y=2))
v22=log(PXY_cond(5)/PXY_cond(6)); %v2(2)=log(P(Xchosen=2|Y=2)/P(Xchosen=0|Y=2))
v21=log(PXY_cond(2)/PXY_cond(3)); %v2(1)=log(P(Xchosen=2|Y=1)/P(Xchosen=0|Y=1))          

U_11_CS=u11;
U_12_CS=u21;
U_21_CS=u12;
U_22_CS=u22;
save('U_CS.mat', 'U_11_CS', 'U_12_CS', 'U_21_CS', 'U_22_CS')

V_11_CS=v11;
V_12_CS=v12;
V_21_CS=v21;
V_22_CS=v22;
save('V_CS.mat', 'V_11_CS', 'V_12_CS', 'V_21_CS', 'V_22_CS')



